##function_call.r

line.color = 'black' ###line color on graphs

##cluster function
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  uj = uj[complete.cases(uj),]  ##modify this function so that it drops NA's become cities don't show up
  
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  return(vcovCL) }

##modified function specifically for polr calls in obs_data_2009.r
 cl.polr   <- function(dat,fm, cluster){
   require(sandwich, quietly = TRUE)
   require(lmtest, quietly = TRUE)
   M <- length(unique(cluster))
   N <- length(cluster)
   K <- 10 #changed it, since polr doesn't generate a variable for rank
   dfc <- (M/(M-1))*((N-1)/(N-K))
   uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
   uj = uj[complete.cases(uj),]  ##modify this function so that it drops NA's become cities don't show up
   vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
   return(vcovCL) }



###cluster function, corrected for small units 
##from http://gmusto.site90.net/bias-reduced-linearization-estimator-in-r/
brl<-function(formula, dataframe, clusterName, method="brl"){
  require(lmtest, quietly = TRUE)
  require(Matrix, quietly = TRUE)
  
  ## standard OLS regression. x=TRUE tells the 
  ## lm function to return the model matrix. 
  ols <- lm(formula,dataframe,x=TRUE)
  
  ## create augmented formula with the cluster variable
  ## and use the new formula to retrieve the cluster
  ## variable from the model.matrix 
  form.clust <- as.formula(paste(deparse(formula), clusterName, sep="+"))
  cluster <- as.matrix(model.matrix(form.clust,dataframe)[,clusterName])
  
  ## store the ols residuals in the vector
  ei <- as.matrix(ols$residuals) 
  
  ## Create the X matrix for the specific regression             
  x <- ols$x 
  
  ## Compute the bread matrix (X'X)^{-1} in our sandwich estimator.                         
  xxm1 <- solve(t(x) %*% x)   
  
  ## Store the number of clusters                
  M <- length(unique(cluster))
  
  ## Store the number of observations
  N <- length(cluster)
  
  ## Store the rank of the X matrix
  K <- ols$rank
  
  # Store the correction for the JK estimator
  dfc <- ((M-1)/M)
  
  "%^%" <- function(x, n){ 
    with(eigen(x), vectors %*% (values^n * t(vectors)))
  }
  
  ## create an empty squared matrix of size ncol(x)
  ## i.e. the number of variables of the model 
  meat=matrix(0,ncol(x),ncol(x))
  
  ## loop over the M clusters
  for (i in 1:M){
    
    ## create an index for the observations in cluster i  
    clusterIndex <- which(cluster==unique(cluster)[i])
    
    ## save the X matrix with observations in cluster i
    Xg <- x[clusterIndex,]
    
    ## generate an identity matrix of size nrow(Xg)
    Ig <- diag(nrow(Xg))
    
    ## compute the matrix Hgg
    Hgg <- Xg %*% xxm1 %*% t(Xg)
    
    ## generate a warning if (Ig-Hgg) is not full rank
    alert<-rankMatrix(Ig-Hgg)
    if (alert < nrow(Xg)){
      warning("no full rank matrix for cluster ", unique(cluster)[i] )
    }
    
    ## compute the residuals correction depending on the  
    ## chosen method
    if (method=="brl") Ag <- (Ig-Hgg) %^% (-0.5)
    if (method=="jk") Ag <- solve(Ig-Hgg) 
    
    ## transform the residual with either the BRL or the JK correction
    rg <- Ag %*% ei[clusterIndex,]
    
    ## compute the meat matrix
    meat <- (t(Xg) %*% rg %*% t(rg) %*% Xg) + meat
  }
  
  ## depending on the method computes the var-cov matrix
  if (method=="brl") vcovBRL <- xxm1 %*% meat %*% xxm1
  if (method=="jk") vcovBRL <- dfc * (xxm1 %*% meat %*% xxm1)    
  
  ## perform t-tests using the var-cov matrix we just computed     
  coef<-coeftest(ols, vcovBRL)
  }


normalize =function (v){
  v/sqrt(sum(v^2, na.rm = T))
  }

##function for converting to upper case
capwords <- function(s, strict = FALSE) {
  cap <- function(s) paste(toupper(substring(s, 1, 1)),
  {s <- substring(s, 2); if(strict) tolower(s) else s},
  sep = "", collapse = " " )
  sapply(strsplit(s, split = " "), cap, USE.NAMES = !is.null(names(s)))
  }
